Voltar ao Artigo
Temperetura RIDE
Descarregar código fonte

Temperetura RIDE

Autor
Afiliação

Carolina Musso

Universidade de Brasília

Data de Publicação

17 de julho de 2025

Resumo

In September 2021, a significant jump in seismic activity on the island of La Palma (Canary Islands, Spain) signaled the start of a volcanic crisis that still continues at the time of writing. Earthquake data is continually collected and published by the Instituto Geográphico Nacional (IGN). …

Palavras-chave

Krigging, Temperatura

Introdução

A temperatura doplaneta lalala

In [1]:
rm(list=ls())
if(!require("pacman")) install.packages("pacman")
Loading required package: pacman
pacman::p_load(
  tidyverse, 
  sf, 
  rio,
  rnaturalearth, 
  rnaturalearthdata, 
  sp, 
  gstat, 
  mapview,
  ggspatial,
  modisfast,
  automap,
  stars
)
In [2]:
## ler os shapefiles com pacote sp
ride_sf <- sf::read_sf("RIDE/.")

ride_sf |> ggplot()+
  geom_sf() +
  scale_fill_viridis_c() +
  theme_minimal() +
  labs(
    title = "Temperatura RIDE",
    subtitle = "Temperatura média do planeta",
    caption = "Fonte: RIDE"
  ) +
  theme(legend.position = "bottom") +
  annotation_scale(location = "br", width_hint = 0.5) +
  annotation_north_arrow(location = "tl", style = north_arrow_fancy_orienteering)

In [3]:
# username <- ""
# password <- ""
# mf_login(credentials = c(username, password))
# 
# 
# roi <- st_as_sf(
#   data.frame(id = "ride_df",
#              geom =  "POLYGON ((-49.36 -17.36, -49.36 -13.10, -45.56 -13.10, -45.56 -17.36, -49.36 -17.36))"),
#   wkt = "geom", crs = 4326
# )
# 
# 
# collection <- "MOD11A1.061" #10:30
# variables <- c("LST_Day_1km")
# time_range <- as.Date(c("2020-10-22", "2020-10-22"))
# 
# urls <- mf_get_url(collection = collection, variables = variables,
#                    roi = roi, time_range = time_range)
In [4]:
# res <-  mf_download_data(urls)
# 
# 
# rast <- mf_import_data(
#   path = unique(dirname(res$destfile))[1],
#   collection = collection,
#   proj_epsg = 4326)
# 
# df <- as.data.frame(rast, xy = TRUE)
# df$lst_c <- df$LST_Day_1km - 273.15
# 
# # 7) Exportar como CSV
# write.csv(df[, c("x","y","lst_c")], "temperatura_geral.csv", row.names = FALSE)
# 
In [5]:
df <- read_csv("temperatura_geral.csv")
Rows: 110038 Columns: 3
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
dbl (3): x, y, lst_c

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
dados_latlong_sf <- st_as_sf(df, coords = c("x", "y"), crs = 4326)

# 2. Reprojetar para o mesmo CRS do município
dados_latlong_sf_proj <- st_transform(dados_latlong_sf, crs = st_crs(ride_sf))

# 3. Filtrar apenas os pontos *dentro* do município
dados_ride_latlong <- st_filter(dados_latlong_sf_proj, ride_sf)  # ✅ correto!
dados_ride_latlong <- dados_ride_latlong[seq(1, nrow(dados_ride_latlong), by = 18), ] #15 x18 16 

# 6. Criar banco com latitude/longitude (graus decimais)
dados_latlong <- dados_ride_latlong %>%
  mutate(
    x = st_coordinates(geometry)[,1],
    y = st_coordinates(geometry)[,2]
  ) %>%
  select(lst_c, x, y) %>%
  as.data.frame() %>%
  select(-geometry) %>% 
  rename(avg_z=lst_c)

# 7. Transformar para UTM (zona 23S - SIRGAS 2000 / EPSG:31983)
dados_utm_sf <- st_transform(dados_ride_latlong, crs = 31983)

# 8. Criar banco com coordenadas UTM (em metros)
dados_utm <- dados_utm_sf %>%
  mutate(
    x = st_coordinates(geometry)[,1],
    y = st_coordinates(geometry)[,2]
  ) %>%
  select(lst_c, x, y) %>%
  as.data.frame() %>%
  select(-geometry)%>% 
  rename(avg_z=lst_c)



# 9. Exportar os dois arquivos
write.csv(dados_latlong, "RIDE_temp_latlong.csv", row.names = FALSE)
write.csv(dados_utm,     "RIDE_temp_utm.csv",     row.names = FALSE)

#set.seed(181326) 
#set.seed(181019)# para reprodutibilidade
set.seed(020279)
dados_utm_sf_amostra <- dados_utm_sf |> 
  rename(temp = lst_c) |>
  sample_n(150)
In [6]:
ggplot() +
  geom_sf(data = ride_sf, fill = NA, color = "black", linewidth = 0.1) +  # polígono em preto
  geom_sf(data = dados_utm_sf_amostra, aes(color = temp), size = 2) +  # pontos coloridos
  scale_color_viridis_c(name = "Temperatura (°C)") +
  theme_minimal() +
  labs(
    title = "Temperatura na região da RIDE",
    subtitle = "Valores extraídos de MODIS (LST)",
    caption = "Fonte: MODIS via modisfast"
  )

In [7]:
v <- variogram(temp ~ 1, data = dados_utm_sf_amostra)
v$dist <- v$dist / 1000  # converte distância para km

fv <- fit.variogram(v, model = vgm(psill = 15, model = "Sph", range = 200, nugget = 3))
plot(v, fv)

fv <- fit.variogram(v, model = vgm(psill = 15, model = "Gau", range = 200, nugget = 3))
plot(v, fv)

fv <- fit.variogram(v, model = vgm(psill = 10, model = "Exp", range = 200, nugget = 2))
plot(v, fv)

In [8]:
# Transformar ambos para UTM zona 23S
ride_utm <- st_transform(ride_sf, 31983)


# ✅ Verificação da variável de temperatura
names(dados_utm_sf_amostra)
[1] "temp"     "geometry"
summary(dados_utm_sf_amostra$temp)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  27.04   31.72   33.41   33.99   35.65   44.91 
# 📈 Variograma empírico e ajuste
vc <- variogram(temp ~ 1, dados_utm_sf_amostra)
vinitial <- vgm(psill = 15, model = "Exp", range = 150000, nugget = 3)
fv <- fit.variogram(vc, model = vinitial)
plot(vc, fv)

# 🔳 Criar grade regular de pontos dentro da RIDE
grid <- st_make_grid(ride_utm, cellsize = 1000, what = "centers")
grid_sf <- st_sf(geometry = grid, crs = st_crs(ride_utm))
grid_sf <- grid_sf[ride_utm, ]

# 🤖 Krigagem
k <- gstat(formula = temp ~ 1, data = dados_utm_sf_amostra, model = fv)
kpred <- predict(k, newdata = grid_sf)
[using ordinary kriging]
kpred_sf <- st_as_sf(kpred)

# 🌍 Projeção para latitude/longitude
ride_latlon <- st_transform(ride_utm, 4326)

kpred_latlon <- st_transform(kpred_sf, 4326)

# 🎨 Plot final com paleta magma + nomes dos municípios\


ggplot() +
  geom_sf(data = kpred_latlon, aes(color = var1.pred), size = 0.6) +
  geom_sf(data = ride_latlon, fill = NA, color = "black",linewidth = 0.1) +
  geom_sf_text(data = ride_latlon, aes(label = NM_MUNICIP), size = 2, color = "black")+   # 
  scale_color_viridis_c(name = "Temperatura prevista (°C)", option = "magma") +
  theme_minimal() +
  labs(
    title = "Krigagem da Temperatura MODIS",
    subtitle = "Interpolação espacial para a RIDE-DF",
    caption = "Fonte: MODIS via modisfast"
  )
Warning in st_point_on_surface.sfc(sf::st_zm(x)): st_point_on_surface may not
give correct results for longitude/latitude data

In [9]:
# 📦 Pacotes
if (!require("pacman")) install.packages("pacman")
pacman::p_load(sf, gstat, stars, leaflet, viridis)

# 🌍 Reprojetar ambos para WGS84 (latitude/longitude)
ride_ll <- st_transform(ride_utm, 4326)
kpred_ll <- st_transform(st_as_sf(kpred), 4326)

# 🎨 Paleta de cores (magma)
pal <- colorNumeric(palette = "magma", domain = kpred_ll$var1.pred)

# 🗺️ Mapa interativo com limites e nomes dos municípios
leaflet() %>%
  addProviderTiles(providers$CartoDB.Positron) %>%
  addCircleMarkers(
    data = kpred_ll,
    radius = 3,
    fillColor = ~pal(var1.pred),
    fillOpacity = 0.7,
    stroke = FALSE,
    label = ~paste0("Temp: ", round(var1.pred, 1), " °C")
  ) %>%
  # Polígono com os municípios (com rótulo interativo)
  addPolygons(
    data = ride_ll,
    color = "white",
    weight = 1,
    fillOpacity = 0,
    label = ~NM_MUNICIP  # altere aqui caso o nome esteja em outra coluna
  ) %>%
  
  # Limite da RIDE destacado em preto
  addPolygons(
    data = st_union(ride_ll),  # união dos polígonos para contorno da RIDE
    color = "black",
    weight = 2,
    fill = FALSE,
    label = "Limite da RIDE"
  ) %>%
  
  # Pontos de krigagem coloridos
  
  
  # Legenda
  addLegend(
    position = "bottomright",
    pal = pal,
    values = kpred_ll$var1.pred,
    title = "Temperatura prevista (°C)",
    opacity = 1
  )

Data & Methods

Conclusion

References